Off-Diagonal Long-Range Order in Solid 4 He 
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Measurements of the moment of inertia by Kim and Chan have found that solid 4 He acts like a 
supersolid at low temperatures. To understand the order in solid 4 He, we have used Path Integral 
Monte Carlo to calculate the off-diagonal long range order (ODLRO) (equivalent to Bose-Einstein 
condensation, BEC). We do not find ODLRO in a perfect hep crystal at the melting density. We 
discuss these results in relation to proposed quantum solid trial functions. We conclude that the 
solid 4 He wave function has correlations which suppress both vacancies and BEC. 
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Whether bulk solid has a region of its equilibrium 
phase diagram with supersolid behavior has been of the- 
oretical interest P, 0, E| for many years. Recently, Kim 
and ChanQ have observed behavior similar to that of 
non-classical rotational inertia (NCRI) at temperatures 
below 0.2K. Analogous to behavior in a superfluid, NCRI 
would be expected in a supersolid, a state of matter in 
which crystallinity and superfluid behavior coexist. Two 
key theoretical quantities in establishing the order of a 
bosonic system are the superfluid fraction (estimated|| 
to vanish) and off-diagonal long range order (ODLRO) 
previously estimated only by variational methods. 

Off-diagonal order (equivalent to a non-zero conden- 
sate fraction) tells us whether atoms at one end of the 
solid are in phase with atoms at the other end of the solid 
and would supply a mechanism for NCRI: for example, 
superfluid 4 He in three dimensions has both ODLRO and 
NCRI. ODLRO 0, is a property of the single particle 
density matrix: 

n(r,r';P) = [ dr 2 -dr n e~ l3E - x (1) 

* Q (r,r 2 ,...,r„)** (r', r 2 , . . . , r„) 

where ^ a and E a are the a th eigenfunction/ value of the 
many body Hamiltonian, Z the partition function, V the 
volume and (5 — l/fc&T. A system is said to be bose- 
condensed if an eigenvalue of the matrix n(r,r') is pro- 
portional to the number of atoms in the system. In this 
paper, we consider a translationally invariant system in 
periodic boundary conditions. This implies that n(r, r') 
only depends on Ar = r — r' and its eigenfunctions are 
momentum states. Considering the occupation of the 
k = momentum state, the condensate fraction is 

no = y5 J drdr ' n ( r > r ')- ( 2 ) 

A system has ODLRO if the condensate fraction is 
greater than zero in the thermodynamic limit. That is: 
n = lim|A r |-»oo n (Ar) > 0. 

ReattojSj showed that a pair-product (Jastrow) trial 
wave function (JWF), commonly used for liquid 4 He, 



has ODLRO. This trial function has the form (R) = 
exp(— 53j<j u ( r ij)) where u (r^) is fixed by minimizing 
the variational estimate of the ground state energy. The 
Jastrow trial function is equivalent to classical system 
with a pair potential v(r): the Jastrow factor maps into 
the classical pair potential (v(r) = 2A:sTu(r)). The phys- 
ical argument for BEC is straightforward: to calculate 
rt(Ar) one displaces an arbitrary atom by Ar; if the av- 
erage amplitude for this is non-zero, the system has BEC. 
The condensate fraction will be non-zero as long as u{r) 
is short-ranged, except for special systems such as hard 
spheres at close packing density. It will also be true if the 
wave function has short-ranged three or higher body cor- 
relation factors. Chester then pointed out that such 
a trial function can be a crystal, since it has the same 
distribution as a classical solid interacting with a pair 
potential. This argument can be used to prove |9J that 
there exist quantum systems with both crystalline or- 
der and BEC. Chester further conjectures that under the 
same conditions, a crystal will have BEC only if there 
exist ground state vacancies (or interstitials). This argu- 
ment has been recently elaborated, (loj 

The Jastrow form of trial function turns out to give 
a seriously incorrect density for freezing ^lj an d other 
properties. As is well-known, melting of classical sys- 
tems occurs when the r.m.s. vibration around a lattice 
site equals about 0.14 of the nearest neighbor distance: 
Lindcmann's melting criterion. However, solid 4 He does 
not melt until much later |l2j. with a Lindemann's ratio 
of 0.30. Assuming Lindemann's criteria for the melting of 
a classical system, no Jastrow trial function can describe 
a solid 4 He near melting, nor can it describe a quantum 
crystal without an intrinsic population of point defects. 

A much better description of solid 4 He is obtained by 
multiplying the Jastrow function by a localization term; 
the Jastrow-Nosanow form|l3( (JN) or "insulator" form: 



p i=i 

Here fa (r) (a Wannier function) is localized about lattice 



2 



site i, and P a permutation of atoms to lattice sites. The 
sum over P projects out the bosonic component of the 
wave function; the projected wavefunction is a "perma- 
nent". In fact, the symmetrization lowers the variational 
energy by an amount prop ortional to the exchange en- 
ergy, calculated in ref. |23 to be about 3 /iK/atom; neg- 
ligible for most properties but crucial for properties such 
as BEC or NCRI. The JN wave function is a very good 
description of the ground state judging by the computed 
energies, Debye- Waller factors, pair correlation function 
and estimated melting density However, the solid or- 
der is put in "by hand" rather than coming about sponta- 
neously. Broken translational symmetry is the hallmark 
of the crystal state[lfij|lfi|. Upon freezing a density wave 
arises which defines the lattice sites, a mean field poten- 
tial and hence <pi(r). It has been shown [lif that this type 
of trial function does not have BEC, assuming that the 
sum of the overlap of <f> with those on other sites is less 
than unity as is expected to be the case. Because of the 
localized functions, the Reatto-Chester theorem on BEC 
does not apply. 

An alternative trial function [ItI Hsl (we will denote 
as a metal, ^ m) is obtained by making the single-body 
function <f>(r) independent of the lattice site, but still 
having lattice symmetry. For electrons, </>(r) would be 
a Bloch function of band theory and would be obtained 
as a solution of a mean field with lattice symmetry. It 
is expected that vacancy- interstitial (VI) fluctuations 
play an important role in supersolidity. In the JN func- 
tion, VI pairs are bound, just as electron-hole pairs are in 
an electronic insulator. In the metallic function, vacan- 
cies will be locally attracted to interstitials but they are 
not bound as pairs; as a result the \&m has both BEC and 
NCRI. Calculations on a similar quantum solid0| have 
found that \E"m has a higher energy than ^ jn\ it costs 
energy to create unbound VI fluctuations. We note that 
calculations [i^ and experiments!!]] on solid 4 He suggest 
that there are no unbound vacancies or interstitials at low 
temperatures. 

Another trial function, the "shadow" wave function |22T] 
(SWF) has been introduced to allow a translationally in- 
variant trial function to have the correct solid order. In 
this case, a single atom coordinate in the Jastrow trial 
function is replaced by a composite object, a linear poly- 
mer, with the inter- and intra- correlation factors being 
variationally determined. Typically the polymer is in fact 
a dimer, with the two ends separated by a distance on 
the order of the inter-atomic spacing. One of the ends 
(the shadow) is integrated over to get the trial function. 
The great advantage of the SWF is that the crystal order 
is spontaneously generated, and one gets accurate prop- 
erties after optimizing the wave function parameters. As 
more monomers are added to the trial function, the func- 
tional form is equivalent to the Feynman-Kacs path in- 
tegral expression for the ground state wave function [2 3|. 
Though it has not been analytically shown, simple argu- 



ments involving displacing the polymer, make it is plau- 
sible that the shadow trial function will have always have 
BEC, as long as the inter-atomic correlation factors are 
short-ranged. This has been verified in recent numerical 
calculations by Galli et al.|24| 

Thus, we have a dilemma, not uncommon with argu- 
ments based on variational wavefunctions: one can have 
several satisfactory trial functions, all of them capable of 
good descriptions of solid helium, but some are BEC and 
some are not BEC. All of the usual properties one uses 
to test the quality of the trial function, e. g. the energy, 
the Debye- Waller factor, the pair correlation function, 
etc. are diagonal in coordinate space (expectations of 
\^>(R)\ 2 ); hence they are unreliable measures of how ac- 
curate the off-diagonal elements needed in Eq. (1) are. 
As we have mentioned above, it is not possible to state 
a general theorem covering whether all quantum crys- 
tals must or cannot have BEC, since the Jastrow trial 
functions can be used to define a Hamiltonian which is 
supersolid. We need a more direct, reliable method to de- 
cide whether the ground state of solid 4 He is BEC and, 
hence, which class of trial function, metal or insulator, 
is appropriate. Calculating n(Ar) in solid bulk 4 He al- 
lows us to establish how good these trial functions are 
for computing BEC in quantum crystals. 

To provide an unbiased answer, we use path integral 
Monte Carlo (PIMC)[2j|, a numerical method that cal- 
culates integrals over the many body density matrix. It 
is ideally suited for this calculation since it can be done 
at finite temperature (under conditions where an experi- 
mental signature of NCRI has been seen), is, in principle, 
exact, and has been validated on many properties|23 of 
liquid and solid 4 He. Most importantly, it is independent 
of a trial wave function bias or any assumption of lat- 
tice. Only the He-He potential enters: a semi-empirical 
form|26| is known to be accurate; in any case, experimen- 
tal results suggest that supersolid behavior is a robust 
phenomenon insensitive to fine details of the interaction. 

In PIMC, based on Feynman's description of super- 
fluid helium, one maps the quantum system onto a clas- 
sical system of "polymers" which can permute onto one 
another. The value of n (| Ar|) is mapped to the end-to- 
end distribution of an open "polymer" in a sea of closed 
"polymers" representing the other helium atoms. If the 
system has ODLRO, the two ends will become separated 
from each other; the condensate fraction is the value of 
n(Ar) at large r. On the other hand, if there is no 
ODLRO, the two ends will remain localized relative to 
each other and n(Ar) will approach at large |Ar|. Su- 
perfluidity (i.e. NCRI), on the other hand, is represented 
in PIMC as long cyclic permutation of paths, for exam- 
ple, as a complete row of atoms shifting down one ele- 
ment. One can immediately see the difference between 
ODLRO and NCRI in a PIMC solid calculation. Sup- 
pose we have N atoms and N lattice sites, so each site is, 
on the average, singly occupied. To perform the ODLRO 
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calculation, one of the atoms paths is opened up, so now 
there will be N — 1 closed polymers and 2 open ends; 
clearly one of the N sites must be doubly occupied if 
we count the open ends as distinct atoms. For NCRI one 
can simultaneously shift all atoms to their new permuted 
sites; there need not be any double occupation, but such 
a simultaneous shift is unlikely for many atoms. 

Techniques^ developed and testedji^ on superfluid 
4 He are used to calculate n(Ar) efficiently, although new 
moves are used to ensure ergodicity in the solid. If one 
simply cuts open one polymer, one rarely finds large val- 
ues of Ar. To guarantee sufficient statistics on the end- 
end distance, we multiply the path integral density by a 
series of importance sampling functions having the effect 
of forcing apart the ends. A series of importance func- 
tions are used to get good statistics for a specific range 
of magnitudes of |Ar| The separate end-to-end distribu- 
tions are then normalized with respect to each other by 
minimizing the x 2 difference between them and using the 
definition ri(0) = I. We perform calculations of n(Ar) 
both by performing an angular average over Ar as well 
as fixing Ar in the basal plane, specifically in the x direc- 
tion. Calculations of ro(Ar) can much more easily access 
the full path space than calculations of the superfluid 
density because the ends of the open polymer can move 
more freely. 

Calculations have been performed in a nearly cubic box 
of 180 particles as well as a series of longer rectangular 
boxes extended in the basal plane ranging from 48 to 144 
atoms. For the longer box, we altered our importance 
function to force the open ends apart only in the x direc- 
tion (in the nearest neighbor direction of the basal plane). 
We examine values of n(Ar) for a range of temperatures 
from Q.IK to 2K in a (periodic) rectangular box of size 
~ 18.5 A corresponding to a density, 0.0286 A~ 3 near the 
experimental melting pressure of 25.3 bars. The random 
walk is started with the atoms in a perfect hep crys- 
talline lattice but no constraints are placed on how the 
atoms can arrange themselves. Since the box is chosen 
to be commensurate with a perfect hep crystal, vacancies 
and interstitials are expected not to exist in the system 
though they can form by a fluctuation. We have verified 
that the results are independent of temperatures below 
2.0K. 

The results of these calculations for T = 0.5K are 
shown in figure 1. Note the curvature at small Ar is de- 
termined by the kinetic energy estimated at 24.0 K/atom. 
For both the spherical average value and the value in the 
basal plane n (|Ar|) decays exponentially for |Ar| > 3A: 
we find no indication of ODLRO (BEC) in bulk solid 4 Hc. 
In fact the rate of decrease for n(x) 1 the SPDM in the x- 
direction, is in agreement with the frequencies for straight 
line winding exchange 5j, shown as the solid line. We find 
oscillations in the computed n(x) reflecting the underly- 
ing crystal lattice. The angular averaged n(r) shows less 
structure than that in the x-direction. Though it is pos- 



sible that n{x) will plateau for a longer distance scale, we 
see no indication of this in the results, nor can we think 
of a physical mechanism that would be responsible |2flj. 
A system with ODLRO will be manifested when one end 
of the open polymer loses knowledge about the other end 
of the open polymer. But we find that pulling apart the 
two ends takes "work" per unit length, no matter how far 
apart the two ends are. This happens because additional 
atoms are displaced from their lattice sites. 

The results of our PIMC calculations are in quanti- 
tative agreement with variational Monte Carlo calcula- 
tions using a shadow wave function [24| for distances less 
than 9 A; they conclude that BEC does exist in the 
ground state of 4 He with a very small condensate frac- 
tion, no = 5 x 10~ 6 at the melting density. However, we 
find that n(Ar) continues to decreases beyond 9 A. 

We believe the reason for the discrepancy with PIMC 
is that the shadow wave function is a technique that 
improves upon the Jastrow wave function, in fact it 
is remarkably accurate up to the second neighbor dis- 
tance. Nonetheless the wavefunction is still built from 
short range correlations. This leads inevitably to a non- 
zero superfluid fraction. We also expect the SWF will 
give a non-zero but very small vacancy concentration pjof . 
PIMC does not make a variational ansatz, with the sole 
assumption that the results at low temperature smoothly 
approach the ground state values. 

Among the proposed quantum solid wave functions, 
the symmetry-breaking Nosanow-Jastrow wavefunction, 
gives results in agreement with PIMC, validating the "in- 
sulator" representation of solid helium. (Further details 
about this picture are in ref. (l^V The implication of 
our calculation and the Reatto-Chester theorem is that 
a quantum solid trial function must include long-range 
correlations. This is most simply done by putting in lo- 
calized functions about lattice sites, reflecting the broken 
symmetry. We have only done PIMC calculations at the 
melting density but we do not expect different behav- 
ior at higher density, since difficulty of exchange grows 
rapidly with density. The result reported here, together 
with the finding of zero superfluid density reported in ref. 
suggest that the mechanism for the measurements of 
Kim and ChanQ involves more than equilibrium prop- 
erties of a perfect 4 He crystal. 

(After the calculation reported here was completed, 
we have learned of a similar PIMC calculation [32| . 
The overall conclusions are similar, however their Monte 
Carlo algorithm and results for n(Ar) differ in important 
ways.) 

This work was supported by NSF under grant no. 
DMR-03 25939 ITR and the fundamental physics pro- 
gram at NASA (NAG-8-1760). Computer time has been 
provided by NCSA and the Materials Computation Cen- 
ter and the Frederick Seitz Materials Research Labora- 
tory (U.S. DOE DEFG02-91ER45439) at the University 
of Illinois Urbana-Champaign. We acknowledge useful 



4 



10' 



10 



10 



10 





1 1 1 1 


1 1 1 


i I i i 


1 




\ 






- 




v\ 




















\ J 
































































V — ■ " 










1 


v . 



10 



FIG. 1: (color online) The single particle density matrix as a 
function of distance (in A), estimated with PIMC, in hep 4 He 
at a density 286nm -3 and a temperature of 0.5K. The red 
curve is the spherically averaged n(r) while the blue curve is 
n(x). The upper (green) curve is a variational MC calculation 
using the shadow wave function ref. |2^ ]. The straight line 
(black) has a slope determined from the exchange calculations 
of ref. |3 • 
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